Astronomy & Astrophysics manuscript no. genbeta-astroph 


© ESO 2008 


February 5, 2008 





Dynamical models with a general anisotropy profile 

M. Baes and E. Van Hese 

Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium 
Received / Accepted 

ABSTRACT 

Aims. Both numerical simulations and observational evidence indicate that the outer regions of galaxies and dark matter haloes are 
typically mildly to significantly radially anisotropic. The inner regions can be significantly non-isotropic, depending on the dynamical 
formation and evolution processes. In an attempt to break the lack of simple dynamical models that can reproduce this behaviour, we 
explore a technique to construct dynamical models with an arbitrary density and an arbitrary anisotropy profile. 
Methods. We outline a general construction method and propose a more practical approach based on a parameterized anisotropy 
profile. This approach consists of fitting the density of the model with a set of dynamical components, each of which have the same 
anisotropy profile. Using this approach we avoid the delicate fine-tuning difficulties other fitting techniques typically encounter when 
constructing radially anisotropic models. 

Results. We present a model anisotropy profile that generalizes the Osipkov-Merritt profile, and that can represent any smooth 
monotonic anisotropy profile. Based on this model anisotropy profile, we construct a very general seven-parameter set of dynamical 
components for which the most important dynamical properties can be calculated analytically. We use the results to look for simple 
one-component dynamical models that generate simple potential-density pairs while still supporting a flexible anisotropy profile. We 
present families of Plummer and Hernquist models in which the anisotropy at small and large radii can be chosen as free parameters. 
We also generalize these two families to a three-parameter family that self-consistently generates the set of Veltmann potential-density 
pairs. These new analytical models are an important step forward compared to isotropic or Osipkov-Merritt models and can be used 
to generate the initial conditions for realistic simulations of galaxies or dark matter haloes. 
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1. Introduction 

Numerical simulations, in particular iV-body or hydrodynamical simulations, have become a major tool to study the structure, 
dynamics, stability and evolution of stellar systems and dark matter haloes. While at first sight this might seem to imply that 
analytical studies of spherically symmetric dynamical models have become less interesting, actually the opposite conclusion should 
be made. The construction of realistic and simple dynamical models is utterly important, as these models act as a reference frame 
or starting point from which to generate the initial conditions for numerical simulations. In this context, it is important to stress that 
the initial conditions need to be generated from the correct distribution function F(r, v) which completely determines the dynamical 
model. Kazantzidis et al. (2004) have recently demonstrated that the correct details of the velocity distribution should be taken into 
account. Simple ad-hoc Jeans dynamical models, i.e. models in which the kinematics are assumed to be Gaussian distributions with 
velocity dispersions derived from solving the Jeans equation, are not sufficient and can lead to erroneous conclusions. The quest 
for analytical dynamical models, which are simple enough on the one hand and realistic enough on the other hand, is hence still 
important. 

In the past few years, many dynamical models have been proposed that are generated by various spherically symmetric potential- 
density pairs. The early-days models mainly represented systems with a constant density core, such as the Plummer model or the 
isochrone sphere (Plummer 1911; Henon 1960). It now appears that many dynamical systems such as galaxies and dark matter 
haloes contain a central density cusp; also for such models a number of representative potential-density pairs have been constructed 
and distribution function have been derived (Jaffe 1983; Hernquist 1990; Dehnen 1993; Tremaine et al. 1994; Hiotelis 1994; Zhao 
1996; Baes & Dejonghe 2002, 2004; Buyle et al. 2007). 

Although the construction of analytical distribution functions hence has continued to develop, most simple models still show an 
important shortcoming: usually the distribution function can only be derived under strict and unrealistic conditions on the anisotropy 
of the velocity distribution. For many potential-density pairs, analytical distribution functions have only been derived under the 
assumption of isotropy. In this case, the phase-space distribution function that supports any density profile can be found as a 
simple integration using the famous Eddington formula. There are a few popular generalizations of the simple isotropic case that 
lead to dynamical models with an anisotropic velocity dispersion. Distribution functions with a constant anisotropy profile can be 
generated with a formula that is very similar to the Eddington formula (Dejonghe 1986; Baes & Dejonghe 2004; An & Evans 2006a). 
Another special category are the Osipkov-Merritt models, engineered independently by Osipkov (1979) and Merritt (1985). In these 
essentially one-dimensional models, the disnibution function depends on the energy and angular momentum integrals only through 
a linear combination of both. The anisotropy profile for these models is peculiar in the sense that the velocity dispersion tensor 
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is isotropic at small radii and becomes completely radial at large radii. The Osipkov-Merritt models were extended by Cuddeford 
(1991) to models in which the velocity dispersion tensor has an arbitrary anisotropy at small radii and becomes completely radial 
at large radii. The transition formulae for all these generalizations of the Eddington formula are sufficiently simple that analytical 
distribution functions of these kinds can be constructed for many simple potential-density pairs. 

Whereas these models are a step beyond isotropic dynamical models, they still do not correspond to the anisotropy that is 
observed both in numerical simulations and real galaxies. Realistic dynamical systems typically have a tendency towards moderate 
to strong radial anisotropy at large radii. Cosmological simulations generally yield dynamical structures that are far from isotropic. 
Dark matter simulations in a cosmological CDM framework typically result in haloes in which the anisotropy gradually increases 
outward to levels of j3 « 0.5 at the virial radius (Cole & Lacey 1996; Colin et al. 2000; Fukushige & Makino 2001; Diemand et al. 
2004). In general, there appears to be a connection between the logarithmic density slope y{r) = dlnp/dlnr(r) and fi(r), both for 
pure dark matter haloes and for structures containing dark matter and baryons (Hansen & Moore 2006). This connection was argued 
to be a natural consequence of the Jeans equation (Dehnen & McLaughlin 2005) and it was shown numerically to be an attractor 
(Hansen & Stadel 2006). 

In hydrodynamical simulations in which gas and stars are taken into account, galaxies also typically have a significant radial 
anisotropy at large radii (e.g. Onorbe et al. 2007). Observational dynamical evidence appears to support these trends. For example, 
detailed stellar dynamical modelling of a large set of elliptical galaxies by Kronawitter et al. (2000) indicates that these galaxies 
generally contain an anisotropy profile that increases from near isotropy in the central regions towards a significant radial anisotropy 
at larger radii. 

The anisotropy at small radii on the other hand is more subtle and depends largely on the dynamical processes that shape 
the nucleus of the system. In particular, the influence of a supermassive black hole, now believed to be present in nearly all hot 
stellar systems, can have a lasting influence on the central anisotropy. Models that contain a supermassive black hole that grows 
adiabatically by slow accretion of material are characterized by an isotropic to slightly tangential anisotropy profile in the central 

region with /3 0.3 (Quinlan et al. 1995). On the other hand, models with a binary supermassive black hole, formed by spiralling 

in due to dynamical friction, typically have much more outspoken tangential anisotropy in the central regions, with /3 1 (Quinlan 

& Hernquist 1997). This enhanced tangential anisotropy results from the preferential ejection of stars on radial orbits during the 
hardening of the black hole binary. Observationally, the study of the anisotropy profile at very small radii is difficult and hampered by 
both limited resolution effects and strong projection effects. Gebhardt et al. (2003) found that the most massive galaxies are strongly 
biased towards tangential orbits in the innermost regions, whereas the lower-mass galaxies have a range of central anisotropies. 

From this body of evidence it is clear that we need to extend our set of simple isotropic or Osipkov-Merritt dynamical models 
to models in which the anisotropy profile has a more general shape, typically rising from isotropy or moderate tangential anisotropy 
at small radii to significant radial anisotropy at large radii. In the general anisotropic case, the transition formulae between density 
and distribution function are rather cumbersome (although mathematically elegant) and involve Laplace-Mellin integral transforms 
(Dejonghe 1986). The goal of the present paper is to describe a method to construct dynamical models with an arbitrary potential- 
density pair and an arbitrary anisotropy profile. We describe the general approach and the mathematical formulation in Section 2. 
In the next two sections we describe a more practical approach to construct such models: in Section 3 we propose a general 
parameterized anisotropy profile and in Section 4 we construct a general seven-parameter family of dynamical components based on 
this parameterized anisotropy profile for which the most important dynamical properties can be calculated analytically. In Section 5 
we apply the results from the previous sections to construct a family of Plummer and Hernquist models with a smooth anisotropy 
profile with arbitrary values at small and large radii. We then extend these two families to a three-parameter family of generally 
anisotropic models that support the Veltmann set of potential-density pairs. Finally, our results are summarized in Section 6. 



2. Construction of general anisotropic models 

2. 1 . The augmented density formalism 

In spherical symmetry, the distribution function depends only on two isolating integrals: the binding energy <5 and the magnitude of 
the angular momentum L per unit mass, 

£> = ft( r ) ~\^~\ v t> (1) 
L — rvj. (2) 

In these expressions if/(r) is the positive binding potential of the system, connected to the density through Poisson's equation 

1 d C^W) = -4,G P (,). (3) 



r 2 dr \ dr 

and v T is the transverse velocity 

Vt 



When we have an expression for the distribution function, all interesting dynamical properties can be derived. In particular, the 
moments of the distribution function can be calculated as 

H2n,im{r) = 2nM F(&, L) v 2 r " vj m+1 dv r dv T . (5) 
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If we want to construct dynamical models for a given density profile in practice, it is easier to consider an alternative characterization 
of a spherical dynamical model, namely the framework of the augmented densities. An augmented density is a bivariate function 
p(if/, r), which expresses the density as an explicit function of both the radius r and the gravitational potential iff. One can demonstrate 
that each augmented density function p(if/, r) completely determines a distribution function F(S, L) and vice versa. There are several 
transition formulae between these two equivalent formalisms, including an elegant formalism based on combined Laplace-Mellin 
integral transforms. 

One of the advantages of the augmented density formalism is that it is straightforward to construct self-consistent dynamical 
models that correspond to a given density profile. Indeed, we just have to solve Poisson's equation, and construct any augmented 
density p(iff) that satisfies the relation 

p(r)=pty(r),r). (6) 

Obviously, there are always infinitely many augmented density functions that satisfy equation (6) for a given density profile p(r). 
These different dynamical models have different moments of the distribution function, and we often are looking for those particular 
dynamical models which have a certain anisotropy profile, defined as 

a-l(r) 

^ = 1 - t^Ty (7) 

2crj(r) 

It appears that the augmented density formalism of dynamical models is a very useful way to achieve this goal. Indeed, another ad- 
vantage of this formalism is that the moments of the distribution function can be derived from the augmented density by subsequent 
derivations and a single integration, 

yn+n Yin + I) r*P r _ 

WM = -Ff7 — 4 (&-f > r*~ Pw.o W, (8) 

Vtt T{m + n) Jo 

where D™ denotes the rath differentiation with respect to x. In particular, the radial and transverse velocity dispersions can be found 
from the density through the relations 

*?(*r) = «J4 = ^ P/^.W, (9) 



PooM, r ) P(fi> r ) Jo 

0*A*.r) = = ^ fV [tW.rjiW. (10) 

MooW n p(iff, r) Jo 

2.2. Construction of general anisotropic models 

Now consider an augmented density that is a separable function of if/ and r, 

p(ifr,r)=mg(r). (11) 

For such models, the dispersion profiles read 

<x^,r) = -i- | (12) 
/W Jo 

^*^ 1+1 2^)mS! md ^' (13) 

such that 

m=-\^{r). (14) 
2 dlnr 

For models with a separable augmented density, the radial velocity dispersion depends only on the (/^-dependent part of the aug- 
mented density, whereas the anisotropy depends only on the r-dependent part. 

Equation (14) offers us in principle the opportunity to construct dynamical models with an arbitrary density profile p(r) and 
an arbitrary anisotropy profile. First we solve equation (14) for g(r). Next we determine the gravitational potential tfr(r) by solving 
Poisson's equation, we invert this relation as r(ifr) and we set 

g(i/f) = g(rW), (15) 
M)=PW)), (16) 

M) = (17) 

then the augmented density 

m,r)=fW)g{r) (18) 

defines the desired model. 
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2.3. Practical construction 

Although we have defined a way to construct an augmented density for a dynamical model with an arbitrary density and anisotropy 
profile, this is of limited use in practice. Indeed, the quantity we really want is the distribution function for the model. In principle, 
the distribution function for each p(ip, r) can be recovered with the standard Laplace-Mellin integral transform formulae. However, 
the augmented density that results from an arbitrary choice of p(r) and f3{r) will often be too complicated to allow an analytical 
integral transform. As a result, we have to compute the distribution function by numerical means. Unfortunately, the numerical 
inversion of equation (5) to obtain the distribution function from the augmented density is in general numerically unstable. Indeed, 
since the density is an integration of the distribution function over velocity space, the density is generally much smoother than the 
distribution function. The inversion procedure of determining / from p hence has the tricky job of unsmoothing the information 
contained in the smoothed density. It hence comes as no surprise that this operation is numerically unstable. For a more precise 
mathematical demonstration of the unstable character of the inversion formulae, we refer to Dejonghe (1986). 

A more promising strategy is to determine a general parameterized model for /3(r) that leads to a relatively simple function 
g(r). For each function f(ifi), we can subsequently construct a dynamical model p(if/, r) that has the same anisotropy profile. Now 
consider a set of base functions ft (if/) that are sufficiently simple such that the augmented density pc(iff, r) = ftiijj) g{r) can be 
converted analytically to a distribution function F( (£, L). It is easy to see that any linear combination of such dynamical models will 
again define a dynamical model with the same anisotropy profile since 



t L t 



g{r). (19) 



So if we can find a suitable linear combination of the base functions such that 

/w=^ = y>/fW, (20) 

we immediately obtain the required distribution function 

F(S,L)=Y j a ( F i (8,L). (21) 

e 

Finding the best linear combination of the components such that the condition (20) is satisfied can in principle be achieved by 
making a formal series expansion and determine each of the expansion coefficients for each component. In practice however, it is 
more straightforward to approximate the dynamical model by a x 2 -minimization routine. This could for example be achieved by 
means of quadratic programming, in which a boundary condition is set that forces the distribution function to remain positive in 
the entire phase space. Practical realizations of this technique are already in use by various dynamical modellers (Dejonghe 1989; 
Kuijken & Merrifield 1993; Merrifield & Kuijken 1994; Gerhard et al. 1998). 

The remaining steps we still have to take are hence (1) create a general parameterized anisotropy profile /3(r) that is flexible 
enough to represent an array of realistic anisotropy profiles and that still yields a simple g{r) function, (2) provide a suitable set of 
base functions or components fy(ijr) for which, in combination with the g(r) function derived above, the corresponding distribution 
function can be computed analytically. We will look into these two steps in the next two sections. 



3. A parameterized anisotropy profile 

As mentioned before, very few analytical dynamical models are known that have a realistic anisotropy profile. Most are either 
completely isotropic, have a constant anisotropy or have an anisotropy profile according to the Osipkov-Merritt type (isotropic 
in the centre and completely radially anisotropic at large radii). In this paper, we aim at the construction of spherical dynamical 
models with a more general and realistic anisotropy profile. The most realistic approach would involve an anisotropy profile that 
depends explicitly on the density slope, as described in the Introduction. However, this would in general not yield an augmented 
density that allows an analytical distribution function. Instead, we will create dynamical models with an anisotropy profile that 
is an explicit function of the radius r. As long as the anisotropy profiles have enough freedom, this will in practice enable us to 
construct dynamical models in which the anisotropy and density slope are at least qualitatively connected. Our specific goal is to 
create dynamical models in which the anisotropy changes monotonically from an arbitrary value fio at the central regions to another 
arbitrary value at large radii, without a priori limitations on the values of /3q or y6 M . To reach this goal, we need to provide a 
functional form that has the desired properties at small and large radii and that can be inverted to give a reasonable simple form 
for g(r). Finding such a form for f3(r) is not obvious. A straightforward candidate that has the desired asymptotic behaviour is an 
exponential profile 



j3(r) =£00-0600- #,) exp 



(22) 



Solving for equation (14) yields, apart from an arbitrary normalization factor, 



g(r) 



-2A,, 



exp 



-2G8 M -j8 )Eii 



(23) 
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It comes as no surprise that such a model will not result in a simple augmented density that can be inverted to an analytical 
distribution function. In order to find a more suitable candidate function for f3(r) we inspire ourselves on the Osipkov-Merritt 
models, where the anisotropy profile reads 

- 2 (r/r. d ) 2 
1 + Wr a ) 2 



resulting in isotropy in the central regions and complete radial anisotropy at large radii. The resulting solution for g(r) reads 

*M = ( 1 + J ■ (25) 
Cuddeford (1991) generalized the Osipkov-Merritt models by considering models with an anisotropy profile 

1 + (r/V a ) 2 

These models also become completely radial at large radii, but at small radii they can assume an arbitrary anisotropy fio- If we 
introduce this expression into equation (14) we obtain the relatively simple expression 

*)-(-) (l^) - (27) 

The obvious way to generalize these results such that the anisotropy does not have to become completely radial at large radii is to 
consider the profile 

M = *+fi-wf t (28) 

1 + (r/r a ) 1 



which corresponds to 



,-2/?o/ 2\-<ft»-A>) 

g(r) = \-\ 1 + ^ . (29) 



In fact, an even more general form for the anisotropy that yields a very similar g function is 



1 + (r/r a ) 2 



with 6 > 0, which corresponds to 



with 



^ I i + ■ (3D 

A = (32) 
o 

The couple (30)-(31) is a suitable choice for our goals. On the one hand this four-parameter model anisotropy contains enough 
flexibility to represent any monotonically varying anisotropy profile (with increasing or decreasing anisotropy). The parameters 
have straightforward effects: /?o and f} x set the anisotropy at small and large radii respectively, r a determines the typical transition 
radius between the two regimes and 6 sets the sharpness by which this transition takes place. On the other hand, the resulting g(r) 
function is sufficiently simple to lead to analytical dynamical components for suitably chosen f(tfr) functions. 



4. A library of components 

The missing link in our approach is a set of base functions or components with which we can construct a linear combination. We 
need to look for functions fi (\jf) such that the augmented density 

r \-2Aw ps\Pe 



PM,r)=fM\T\ 1 + ' (33) 



gives rise to an analytical distribution functions F[(£, L). We will drop the subscript t in the remainder of this section in order not 
to overload the notations. 
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4.1. Power law components 

The most obvious candidate components are those in which f(ifr) is a power law, 



.26 )Ps 
26 



(34) 



where we limit ourselves to systems with a finite potential well if/Q = iff(0). A finite total mass is obtained if p + 2/3^ > 3. In 
order to calculate the distribution function corresponding to this augmented density, we apply the general Mellin-Laplace inversion 
formulae. After some algebra (see Appendix A), we find that it can be expressed as a Fox H-function: 



F(S,L) 



A) 



M(2n^?V 6T(-p s ) \t/r 



m + p) (sy- y2 u 

H 2,2 



2rl& 



(i-^a.(p-M)i 
(-f.a.(o.i) 



(35) 



As we show in Appendix A, a sufficient (but not necessary) condition to yield a well-defined distribution function, i.e. continuous 
and non-negative, is 6 < 1 or 6 = 1 and p - | + yS 5 > 0. 

For practical purposes, a series expansion is more useful. After some lengthy algebra, one obtains the fairly simple expression 
(see Appendix B) 

1 < < 2 ^ 0+M 



F(G,L) 



p T(l+p) IS 



M(2^(Ao) 3/2 \ifro 



P-3/2 



x < 



z 



^U/F(l -fo + k8)T(p~\ +p Q -kS) \2rl6 



'A 



l 



I- 



-k6 



for L 2 < 2r 2 &, 



for L 2 > 2r 2 8. 



(36) 



T(l -jSc - kS)T(p - \ +p„ + fay) \2r 2 S) 

Apart from the distribution function, most other interesting dynamical properties can be calculated analytically for this set of 
components. Particularly simple are the velocity dispersions, which we find through (12), 



o 2 r ty,r) = 



l+p' 



~2 (l v ~2 (i v (l-^o) + (l-^)(r/r a ) M _ 2 



1 + (r/rj 



(37) 
(38) 



4.2. Extension to cuspy components 

This set of power law components presented in the previous subsection is very adequate to fit a broad class of models. It is however, 
not fit to construct dynamical models for systems that have a central density cusp. For such models, the isotropic distribution function 
diverges in the limit £ — ♦ if/Q (see e.g. Hernquist 1990; Dehnen 1993; Tremaine et al. 1994; Baes et al. 2005). This behaviour cannot 
be represented with the current set of components. Therefore we extend the set of models defined by the augmented density (34) to 

(ifr\ p l FVlrY 7 *! r 2s \ ps 

***-"(*) ('-IJUJ r?) - (39) 

with two additional parameters q < and s > 0. This seven-parameter family of dynamical components contains enough freedom 
in the density profile to form a good family of base functions. If we expand this expression as a power series in iff, 



p(ifr, r) = p 



-2/3„ 



1 + 



J18\P6 



.26 



(40) 



we see that we obtain a series of positive terms that all have the form (34). Thus, our condition q < is sufficient to obtain a 
well-defined distribution function, which we can write down immediately, 



F(8,L) = 



A) 



M(27T(Ao) 3/2 



f(-iyj q ) ni+p+ 



6) Wo) 



P+js-3/2 



H 



2rl& 



(-$.*).<P.D 



(41) 



or explicitly 
F(G,L) = 



A) 



M(2wi/r ) 3/2 



!<-»<(«) 



F(l +p + js) 



p+js-3/2 



s 

k=0 

CO 

Z 

l. k=0 



-Po+kS 



k /T(l -fo+k6)T(p + js-\ +/3 -kd) \2r£6 

'fie\ ! ( JM 

kjT(l-p oa -kS)T(p + js-^+Poo+ks) \2r a 2 £/ 



-Poo-k6 



for L 2 < 2rl&, 
for L 2 > 2rl&. 



(42) 
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Since the summations in j converge for every fixed value of k, this distribution function is indeed well-defined. Other dynamical 
properties, such as the moments of the distribution function, can be derived in the same way. For the velocity dispersions we find 



s Wo 



i/r s 



l+p 



1 + q 



\26 



(l- J 8o) + (l-A :o )(r/r a ) 2 „ 

— cr;(^,r), 



1 + (r/rj 



26 



(43) 
(44) 



with B x {a, b) the incomplete Beta function. 

5. Self-consistent analytical models 

In the previous Section we have described a practical method to construct analytical self-consistent dynamical models that generate 
any spherical density profile and with an anisotropy profile defined by the general parameterized function (30). We are applying 
this technique to construct dynamical models matching the density and anisotropy profiles that have come out of cosmological 
simulations of dark matter haloes, in order to gain insight into their phase-space structure (Van Hese, Baes & Dejonghe 2007). 
In the remainder of this current paper we will focus on simple analytical dynamical models, by investigating whether the global 
seven-parameter family of dynamical models defined by the augmented density (39) contains any simple self-consistent models 
(hence without the need to make a linear combination). Thus we will look for potential-density pairs [p{r), <fr(r)] that satisfy both 
Poisson's equation and the relation 

p(r)=p(if,(r),r), (45) 

with p(i[/, r) the seven-parameter augmented density from equation (39). The existence of such analytical dynamical models would 
be a great step forward in our quest for simple but realistic dynamical models that can be used as a framework in which to initiate 
detailed numerical simulations. 



5.1. A family of anisotropic Plummer models 

One of the most obvious candidates is the Plummer model, as this model has a rather straightforward potential-density pair 

GM 



VP" 



+ r L 



P(r) = 



3M 



2\-5/2 



1 + 



Anb 3 y ' b 2 . 

When we combine these functions with the expressions (39) and (45), we obtain the condition 



3M / + ij_ 



-5/2 



Po 1 + 



b 2 



■p/2 



1-1 



b 2 



s/2 



-2/Ju 



1 + 



.25\-(ft=-A>)/<5 



.26 



A straightforward solution is obviously given by 



3M 

P = 5, 
g = 6 =8< X , = Q, 



which yields the isotropic Plummer model, defined by 



3M 



4nb 3 \GM) 

3 ( 2b& 
In 3 (GMb) 3 ' 2 \GM 



7/2 



(46) 
(47) 

(48) 

(49) 

(50) 
(51) 

(52) 
(53) 



Our goal, however, is to determine the most general subspace of the (p, q, s,Bq,B x , 6, r a ) parameter space such that the condition (48) 
is satisfied for all r. In particular we aim for a subspace of the parameter space that has no restrictions on Bq and B^, such that we 
obtain a family of dynamical models with an arbitrary anisotropy at small and large radii. It is obvious from equation (48) that, in 
order to find non-trivial solutions, we have to set 



s = 25 = 2, 
r« = b. 



(54) 
(55) 
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We then obtain the expression 



3M 



2 \Po-q 



Anb^pQ \b 2 

These expressions are identical for all values of r if 



2\-5/2+p/2+<jr +J S„-/3 



1 + 



b 2 



3M 



p = 5-2p oa , 
q = Po- 



(56) 



(57) 

(58) 
(59) 



We now have constructed a general two-parameter family of self-consistent Plummer models with the augmented density profile 



p(i]/, r) = 



3M 



4nb 3 \GM 



btfr\ 



I bifr \ 2 1^° / r 2 \ _jS ° 

1 ~\gm) [b 2 ) 



2\-08=.-/8d) 



1 + 



b 2 



(60) 



From the condition q < it follows that the central anisotropy has to satisfy /3o < 0. This is in fact a special case of the cusp 
slope-central anisotropy theorem of An & Evans (2006b): a model without a density cusp cannot have a radial velocity anisotropy in 
the centre. The anisotropy at large radii can take any value p^ < 1. By construction, the anisotropy profile of this family of models 
reads 

fob 2 +Poor 2 



(61) 



b 2 + r 2 

which indeed tends towards /3q at small radii and towards j6 M at large radii. The radial velocity dispersion profile can be written as 



GM I r 

~ib\b 2 



2\~Po 



2x 5/2-08„-/3 o ) 
l + V 2 



B_^_{3-p„,\+Po). 



At large radii, the radial velocity dispersion profile shows a r 1 behaviour for all values of the parameters 0q and 



1 GM 

6 - 2^00 r 



(62) 



(63) 



The asymptotic behaviour at small radii is 



0?(r) 



GM 1 



b 1 + /3o b 2 



GM 



_JK3-/fc.,l +A ) p 



2\-Po 



if A) < -1, 



ifjfio > -1. 



(64) 



Except for the models that are isotropic in the centre, the radial velocity dispersions hence always tend to zero at small radii. The 
asymptotic behaviour of the tangential velocity dispersions cr g (r) = (r v (r) follows immediately. At large radii we obtain 



2 M i ( , GM 
°W) = o-Jr) 



whereas at small radii 



oiir) = crl(r) 



( _GM 1 -fa r 2 
V 1 +yS 
1 -Po GM 



5(3 - jffoo. 1 + Po) 



b 2 



apo <-h 



if ft, > -1. 



(65) 



(66) 



In the top panels of Fig. 1 we plot the radial and tangential velocity dispersions for a set of anisotropic Plummer models. The three 
models in this figure all have the same anisotropy p m = j at large radii, but a different central anisotropy. Apart from the radial and 
transversal velocity dispersion, we also plot the projected velocity dispersion cr p (R) on the plane of the sky, defined through 



d2 

1 " -jftr) 



r 

Pp (K)o*(B) = 2 

JR 

Here, p p (R) is the projected mass density on the plane of the sky, 

p(r) rdr M b 2 



p(r) cr 2 (r) r dr 



y/r 2 -R 2 



P P {K)=2 f 

JR 



V r 2 _ R 2 7T (b 2 +R 2 ) 2 ' 



(67) 



(68) 
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Fig. 1. Top. Radial (red), tangential (green) and projected (blue) velocity dispersion profiles for three Plummer models with 
anisotropy values Bq and B^ displayed in the figures. Bottom. The distribution function of these models, represented by isoprobabil- 
ity contours in turning point space. High values are indicated by red contours, low values by yellow contours. In all plots, we have 
used normalized units with G — M = b — I. 



The projected velocity dispersion always reaches a finite value in the centre and has an R dependence at large radii, 

2/ 



3;r 6-5/3 M GM 
? ( ' ~ 128 3-B x R " 



The distribution function for our set of Plummer models can be written as 
F(F n= 3 1 y v //?o\ r(6 + 2;-2/U /fr£\ 



7/2+2,-2,8..,. 



9 L 2 \ 



where the function H(a, b, c, d; x) is a function defined as (Dejonghe 1986) 

1 -b,c 



R(a,b,c,d;x) = G 



2.2 



a, 1 - d 



= H 



2.2 



(l-M),(c,l) 

(a,l),(l -d,l) 



with Gp'o (z) the Meijer G-function. This H-function can conveniently be expressed as 



M.(a,b,c, d;x) 



T(a + b) 



T(c — a) F(a + d) 

T(a + b) (I 
{ T(d-b)T(b + c) [x 



x° iF\(a + b, 1 + a — c; a + d; x) 



iF\ [a + b,\ + b - d;b + c; 



1 



if x < 1, 



if x > 1 . 



(69) 



(70) 



(71) 



(72) 



In the bottom panels of Fig. 1 we plot the distribution function as a contour plot in the turning point space for the same models 
as the upper panels. The change in anisotropy from tangentially anisotropic at small radii to radially anisotropic at large radii can 
easily be seen in the slope of these contours. 

An interesting subset of models in our two-parameter family of Plummer models is the one-parameter family with Bq = 0. These 
models are isotropic in the inner regions and become anisotropic at large radii. Their augmented density is given by 



P(<A, r) 



3M 



4tt£ 3 \GM 



bl_\ 5 



1 + 



b 2 



(73) 
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and in the expression for the distribution function (70) only the term corresponding to j — remains 

eyp r> 3 1 r(6-2 J 6 co ) /fc£\ 7/2 - 2/j - mr / A 9 L 2 \ 

F{S,L)= ,~,„,v»n wo ; — 7^77 H 0,j8cc,t -2y8oo,l; — r- . (74) 



This subfamily of our current set of Plummer models was already presented by Dejonghe (1987). Most of the kinematical properties, 
including the projected properties such as dispersions and higher-order moments of the line profiles, can be calculated completely 
analytically. 



5.2. A family of anisotropic Hernquist models 

Another extremely popular and simple potential-density pair is the Hernquist model, defined by 

, s GM 
if,(r) = 



b + r 
■«r r(b + r) 



3 ' 



(75) 
(76) 



Contrary to the Plummer model, this model has a central r~ l density cusp and a more realistic r~ 4 behaviour at large radii. We can do 
the same experiment for the Hernquist model as we did for the Plummer model. If we combine the potential-density pair (75)-(76) 
with expression (39) and (45), we obtain 



M 

1Kb 1, 



afKf-K)>K)~ 



-200 



r ay\-Gs»-A))/<5 

1 + 



Similarly as for the Plummer model, it is clear that we will only be able to find a general non-trivial solution if we set 

s = 26 = 1, 
r a = b. 



This yields the equation 



from which we find 



M 



i r \-i-?+2#> / r \ 



2nb 3 pQ 



M 



Po 



2nb 3 ' 
p = 4-2p 00 , 
q = 2J3 -l. 

We have now defined a two-parameter family of self-consistent Hernquist models with augmented density 



^ = 2^[^j [ l -cM (b) K) 



v-2(ft»-/%) 



ill) 

(78) 
(79) 

(80) 

(81) 

(82) 
(83) 

(84) 



The parameter can assume all values, whereas the central anisotropy /3q is limited to A < i, in agreement with the cusp slope- 
central anisotropy theorem of An & Evans (2006b). By construction, the anisotropy profile of this family of Hernquist models 
reads 

= A*±A^ 

b + r 

which has the desired behaviour at small and large radii. The radial dispersion profile reads 

GM I r^-Wo I r \3-208„-j8 o ) 



^ (r)= TU ( 1 + b) ^v<5-2A,.2/i„). 
At large radii, the radial velocity dispersion profile falls as 



2, , 1 GM 

cr r {r) ~ - — — , 

5 - 2/3^ r 



(86) 



(87) 



whereas the asymptotic behaviour at small radii is 

GM 1 r 



b 2fio b 
GM 
~b~ 



5(5-2^,2^)^) 



l-2ft 



if A) < 0, 
if A) > 0. 



(88) 
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Fig. 2. Top. Radial (red), tangential (green) and projected (blue) velocity dispersion profiles for three Hernquist models with 
anisotropy values Bq and B^ displayed in the figures. Bottom. The distribution function of these models, represented by isoprobabil- 
ity contours in turning point space. High values are indicated by red contours, low values by yellow contours. In all plots, we have 
used normalized units with G — M = b — I. 



The radial velocity dispersion hence always disappears in the centre, except for the models with the largest allowed central 
anisotropy (Bo = j) where it reaches a finite value. For the asymptotic behaviour of the tangential velocity dispersions crg(r) = cr^r) 
at large radii we obtain 



whereas at small radii 



oi(f) = o-\(r) 



( GM 1 -B r 



1 -Be, GM 
5- 2/3^ r ' 



b 2B b 

GM (r\ l 
(l-Bo)B(5-2B x ,2B )[-j 



v b 



if Bo < 0, 
if B > 0. 



The projected velocity dispersion on the plane of the sky has a R 1 dependence at large radii, 

^ ' ~ \5n 5 - 2fi„ R ' 
The distribution function can most conveniently be written as a series of hypergeometric functions 



(89) 



(90) 



(91) 



F(S,L) 



1 



1 



b8 



5/2-3A. 



(2n) 5 ' 2 (GMbfl 2 \GM) 

" lB s \ T(5-2B 0O ) 



^fe/r(i-A + |)r(|-2^ +A -|)V2fc 2 £ 



2 FA5 - 2/3^,1 -2Bo;--B co + --—\ (92) 
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ifL 2 < 2b 2 8, and as 

F( fi ,L)= 1 ' ' b& 



(2tt) 5 / 2 (GMb) 3 ' 2 \GMj 
X 



\5/2-2ft» 



F(5-2/U / L 2 \^- k/2 ( 7 * fc6\ 

if L 2 > 2b 2 &. These sums only contain a finite number of terms if fig is a positive integer number, i.e. when (J3q - fico) is a positive 
integer or half-integer number. This particular subset of models, in which the outer regions are always more tangentially anisotropic 
than the central regions, has already been discussed by Baes & Dejonghe (2002). 

In a similar manner as for the Plummer model, we plot in Fig. 2 the velocity dispersions and the distribution function for three 
Hernquist models with the same anisotropy p x but different anisotropy (3q. 

5.3. Generalization to a three-parameter family of anisotropic Veltmann models 

It is well-known that the Plummer and the Hernquist potential-density pairs can be generalized to a one-parameter family of models 
characterized by 

GM 

t (r) = T71 ^yvy, (94) 

_ (1 + A)M b A 

~ 4^ r 2-X (UX + r A.\2+\IA ' ^ ^ 

This potential-density pair was first described by Veltmann (1979) and is a special subset (the a-models) of the general set of 
potential-density pairs considered by Zhao (1996). This potential-density pair recently regained much interest because it supports 
dynamical models that are hypervirial, i.e. in which the virial relation is not only satisfied on a global but also on a local level (Evans 
& An 2005; Iguchi et al. 2006; Sota et al. 2006). The parameter A, lying in the range < A < 2, determines the slope of the central 
density cusp. We easily recognize the Plummer model with A = 2 as the only core-density member of the family and the Hernquist 
model as the model with A = 1 . 

We can now repeat the same exercise as for the Plummer and Hernquist models. After a little bit of algebra, we find that the 
parameters 

Po = "W (96) 

p = 3 + A-2/3 0O , (97) 
+ (98) 

A 

s = 26 = A, (99) 
r- d = b, (100) 

are the general solution for the condition of self-consistency. Notice that the initial condition q < implies /?o < 1 - A/2, which is in 
correspondence with the cusp slope-central anisotropy theorem (An & Evans 2006b). In other words, for this family the condition 
q < is also necessary to yield physical models. In this manner we have constructed a three-parameter family of dynamical models 
defined by the augmented density 



HW ' \ntf> \GM) I \GMj J \b A ) \ b A J 

This augmented density self-consistently supports the one-parameter potential-density pair (94)-(95) and has by construction an 
anisotropy profile 

b A + r A 

which varies smoothly from fio in the centre towards fico at large radii. The radial velocity dispersion can be written as 



l+2(8o-l)M. A .-2/3 /A , i\-W«-M/A 



For general values of A, the distribution function cannot be simplified and should be taken as in formula (41) with the values (96)- 
(100). For rational values of A however, distribution function simplifies to a sum of generalized hypergeometric functions. One such 
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case is the model with A — 1/2, for which we obtain the potential-density pair 

GM 

W) = -2, (104) 

3M 

P(r) = s ~ ~ 4 - (105) 

87r r 3 / 2 ( V£ + VF) 

This model has a central density cups with a slope p(r) oc r~ 3//2 , which has been obtained by numerical simulations for dark matter 
haloes in a CDM cosmological model (Moore et al. 1998). 



6. Discussion and conclusions 

In this paper, we have described a method to construct self-consistent spherical dynamical models that have an arbitrary density and 
an arbitrary anisotropy profile. While the general formulation is not so useful since the inversion formulae are rather complicated 
and numerically unstable, we propose an approach with a parameterized anisotropy profile. We put forward a very general four- 
parameter anisotropy profile in which the four parameters control the anisotropy at the centre and at large radii, the typical transition 
radius between the two regimes and the sharpness in which this transition takes place. This parameterized function can therefore 
represent a large number of monotonically varying anisotropy profiles. 

Based on these anisotropy profiles, we present a seven-parameter set of dynamical components that all have the same anisotropy. 
For each of these components the augmented density p(tfr, r) is a sufficiently simple function so that the corresponding distribution 
function can be written analytically (as a double series). For many reasonable choices of the parameters, the distribution function 
can be written as a sum of generalized hypergeometric series. The lower-order moments can be expressed by simple analytical 
functions - for example, the velocity dispersions of all components can be expressed as incomplete Beta functions. 

Once the anisotropy profile has been fixed, the construction of a dynamical model with a given potential-density pair consists 
of a determination of the best linear combination of the various components. This can be achieved with a quadratic programming 
approach, a technique already in use by various dynamical modellers (Dejonghe 1989; Kuijken & Merrifield 1993; Merrifield & 
Kuijken 1994; Gerhard et al. 1998). This dynamical modelling technique is rather advanced: dynamical models can be constructed 
by fitting a variety of data, such as any moments of the distribution function, the projected moments along the line of sight, the 
full line profiles or even complete spectra. One could argue what the novelty of our approach is compared to the general powerful 
techniques already presented by these authors. The use of the quadratic programming technique is very general and does not require 
all components to have the same anisotropy profile. Any possible combination of components can in principle be combined. The 
main strategy usually consists in choosing the components in such a way that the distribution functions are simple enough; one can 
for example choose only Fricke components which are double power laws both in augmented density and distribution function. As 
long as the data base of possible models is large enough, such components allow to fit basically any dynamical model since a linear 
combination of two components with a different density profile and different constant anisotropy results in a model with a non- 
constant anisotropy. However, practice has learned that it is generally difficult to construct models with a strong radial anisotropy 
at large radii. The reason for this problem is that one needs to populate the model with radial orbits that reach large radii. Since 
such orbits also contribute to the density at small radii, it requires a delicate fine-tuning of the different components to both satisfy 
the density and anisotropy constraints at small radii while still retaining the radial anisotropy at large radii. In our present case, this 
problem does not surface since we only have to fit the function f(i/f), or equivalently the density profile. Any combination of the base 
components will automatically have the correct anisotropy profile. This was actually one of the motivations for the construction of 
this set of dynamical components. 

Apart from presenting an approach to construct detailed dynamical models with a given density and anisotropy profile, our study 
has also lead to a set of simple one-component dynamical models. The models we have generated are self-consistent realizations 
of the Plummer and Hernquist models. We have also generalized these models to a complete set of models that self-consistently 
generate the family of Veltmann (or generalized Plummer) potential-density pairs. Contrary to most analytical dynamical models 
that have appeared in the literature, the anisotropy profile of these new dynamical models can be chosen arbitrarily in the sense that 
we can choose the anisotropy at both small and large radii. Both numerical simulations and observational evidence clearly indicates 
that the outer regions of galaxies and dark matter haloes are typically mildly to significantly radially anisotropic and that the inner 
regions can be significantly non-isotropic. Realistic simulations should use this information and not impose unrealistic constraints 
on the anisotropy. We therefore believe that these new analytical models are an important step forward compared to the isotropic or 
Osipkov-Merritt type dynamical models that are often used. We encourage numerical astrophysicists to use these new distribution 
functions to generate the initial conditions for their simulations of galaxies or dark matter haloes. Numerical implementation of the 
most important dynamical properties of these models can be obtained from the authors. 

We are currently employing the modelling technique described in this paper to construct dynamical models for dark matter 
haloes (Van Hese, Baes & Dejonghe 2007). We are using constraints on the density and the anisotropy profile of these haloes from 
detailed A^-body simulations. This will be useful to investigate the phase space structure of dynamical models and will also provide 
the community with an advanced toy model for dark haloes from which more simulations can be initiated. This nicely illustrates 
that the rise of numerical simulations does not erase the need for analytical work, but rather that there is a symbiosis between both 
approaches. 



14 



M. Baes & E. Van Hese: Dynamical models with a general anisotropy profile 



References 



An, J. H., & Evans, N. W. 2006a, AJ, 131, 782 

An, J. H., & Evans, N. W. 2006b, ApJ, 642, 752 

Baes, M., & Dejonghe, H. 2002, A&A, 393, 485 

Baes, M., & Dejonghe, H. 2004, MNRAS, 351, 18 

Baes, M., Dejonghe, H., & Buyle, P. 2005, A&A, 432, 411 

Buyle, P., Hunter, C, & Dejonghe, H. 2007, MNRAS, 375, 773 

Carlberg, R. G., et al. 1997, ApJ, 485, L13 

Cole, S., & Lacey, C. 1996, MNRAS, 281, 716 

Coh'n, P., Klypin, A. A., & Kravtsov, A. V. 2000, ApJ, 539, 561 

Cuddeford, P. 1991, MNRAS, 253, 414 

Cuddeford, P., & Louis, P. 1995, MNRAS, 275, 1017 

Dehnen, W. 1993, MNRAS, 265, 250 

Dehnen, W., & McLaughlin, D. E. 2005, MNRAS, 363, 1057 
Dejonghe, H. 1986, Phys. Rep., 133, 217 
Dejonghe, H. 1987, MNRAS, 224, 13 
Dejonghe, H. 1989, ApJ, 343, 113 

Diemand, J., Moore, B., & Stadel, J. 2004, MNRAS, 352, 535 
Evans, N. W., & An, J. 2005, MNRAS, 360, 492 
Fricke, W. 1952, Astronomische Nachrichten, 280, 193 
Gebhardt, K., et al. 2003, ApJ, 583, 92 

Gerhard, O., Jeske, G., Saglia, R. P., & Bender, R. 1998, MNRAS, 295, 197 

Gradshteyn, I. S. & Ryzhik, I. M. 1965, Table of integrals, series and products (New York: Academic Press, 1965, 4th ed., edited by Geronimus, Yu.V. (4th ed.); 

Tseytlin, M.Yu. (4th ed.)) 
Fukushige, T., & Makino, J. 2001, ApJ, 557, 533 
Hansen, S. H, & Moore, B. 2006, New Astronomy, 1 1, 333 

Hansen, S. H, & Stadel, J. 2006, Journal of Cosmology and Astro-Particle Physics, 5, 14 
Henon, M. 1960, Annales d'Astrophysique, 23, 474 
Henon, M. 1973, A&A, 24, 229 
Hernquist, L. 1990, ApJ, 356, 359 
Hiotelis, N. 1994, A&A, 283, 783 

Iguchi, O., Sota, Y, Nakamichi, A., & Morikawa, M. 2006, Phys. Rev. E, 73, 0461 12 
Kuijken, K., & Merrifield, M. R. 1993, MNRAS, 264, 712 
Jarre, W. 1983, MNRAS, 202, 995 

Kazantzidis, S., Magorrian, J., & Moore, B. 2004, ApJ, 601, 37 

Kronawitter, A., Saglia, R. P., Gerhard, O., & Bender, R. 2000, A&AS, 144, 53 

Louis, P. D. 1993, MNRAS, 261, 283 

Mathai, A. M. 1993, A Handbook of Generalized Special Functions for Statistical and Physical Sciences, Oxford University Press, Oxford. 
Merrifield, M. R., & Kuijken, K. 1994, ApJ, 432, 575 
Merritt, D. 1985, AJ, 90, 1027 
Merritt, D. 1985, MNRAS, 214, 25P 

Moore, B., Governato, F, Quinn, T., Stadel, J., & Lake, G. 1998, ApJ, 499, L5 
Ofiorbe, J., Dorm'nguez-Tenreiro, R., Saiz, A., & Serna, A. 2007, MNRAS, 376, 39 
Osipkov, L. P. 1979, Pis ma Astronomicheskii Zhurnal, 5, 77 
Plummer, H. C. 1911, MNRAS, 71. 460 

Quinlan, G. D., Hernquist, L., & Sigurdsson, S. 1995, ApJ, 440, 554 
Quinlan, G. D., & Hernquist, L. 1997, New Astronomy, 2, 533 

Sota, Y, Iguchi, O., Morikawa, M., & Nakamichi, A. 2006, Progress of Theoretical Physics Supplement, 162, 62 

Tremaine, S., Richstone, D. O, Byun, Y.-L, Dressier, A., Faber, S. M., Grillmair, C, Kormendy, J., & Lauer, T. R. 1994, AJ, 107, 634 

Van Hese, E., Baes, M., & Dejonghe, H. 2007, in preparation 

Veltmann, U. I. K. 1979, AZh, 56, 976 

Zhao, H. 1996, MNRAS, 278, 488 



Appendix A: Derivation of the distribution function 

We now present the detailed calculation of the distribution function corresponding with the augmented density (34). We have to 
consider several cases, depending on the parameter 

A.1. Case 1:/} s is a natural number 

First, we consider the special case where f} s is zero or a natural number. Then the augmented density is simply a finite sum of 
positive Fricke components (Fricke 1952; Henon 1973) 




(A.l) 



which leads immediately to equation (36). 
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A.2. Case 2: /3 S < 

When Ps < we can apply the Laplace-Mellin formalism: the connection between the distribution function and the augmented 
density is then (Dejonghe 1986) 

2*7/2 £(3-i7)/2 



£ M{F(8,L)} 



. - £ M{p(if,,r)}. 
M{2nyi l p n _ 2j ^-tf r->n 



Since the augmented density is a separable function of if/ and r, the transforms can be calculated separately: 



Jo 



r(i + P ) 



and 
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r->!7 



1 + 



26 



-r 

Jo 



r»- 1 



26\-W™-MI« 



1 + 



26 



dr = — B 

25 



77 - 2j8 2/3 M - 77 



where 77 lies in the convergence strip 2/?o < 77 < 2/3 M - Thus 



r a. r*YP m - P^' 1 r " r(1 + PK l v- 2 A> ~ g 
( ),_ M<(2^/22tf r i_2) I 25 ' 25 
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<l-i7)/2-p 
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The inversion of the Laplace transform is straightforward, 
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27T7 



£ -1 \^-n)l2-p\ 
f-tS ^ ' 



e f^(l-„)/2- / , d ^ 
-i<* I ' (/) 



gp-(3-77)/2 



lz2 
2 



leaving us with the inversion of the Mellin transform 

Po 



F(8,L) = 



M(27 n Ao) 3/2 5r(-/3 tf )\ ( A„/ 
po r(i 



Milnxjjofl 2 5Y( 
Taking into account the definition of the Fox //-function, 
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ijinji 
H P,q 
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the distribution function can be written as 

Po 



or equivalently 



F(S,L) 



F(8,L) = 
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(A.2) 



(A.3) 



(A.4) 



(A.5) 



(A.6) 



(A.7) 



(A.8) 



(A.9) 



(A. 10) 



The integration can be performed along three possible paths C, which are equivalent: if the integral converges for more than one of 
these three paths, then the result is the same. If the integral converges for only one path, then that is the only one to be considered. 
The valid contours are 

- C is a path C\ from -700 to 700, such that the poles of Tf^y^J lie on one side and the poles of rf^pj lie on the other side. 
The convergence is absolute if 5 < 1 or if 5 = 1 and p - \ + fig > 0. This condition is found by investigating the more familiar 
criterion for Meijer G-functions (see Appendix B). In addition, if 5 = 1 and p + ~ + /3g > the integral is semi-convergent for 
L 2 + 2r\& (Dejonghe 1986, Appendix). 

- C is a loop C2, starting and ending at -00, that encircles the poles of r( ^y^J once in the positive direction and none of the poles 
of T J. The integral then converges if L 2 < 2r 2 £. 

- C is a loop C3, starting and ending at +00, that encircles the poles of r(y^) once in the negative direction and none of the 
poles of F f^pj. The integral then converges if L 2 > 2r 2 &. 

It can easily be seen that the convergence criterion for path C\ ensures a well-defined distribution function, i.e. continuous and 
non-negative. Indeed, in this case the contour is a line parallel to the imaginary axis from sq - 700 to sq + 700 with /3q < So < On 
this path the real part of the Gamma functions is positive (where the condition p + 2/3 M > 3 ensures that p-^ + s>0by choosing 
i(j sufficiently close to /3 M ). Hence, the distribution function is non-negative everywhere. 
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A3. Case 3:j3 s > and not a natural number 

In the case of fig > and not a natural number, the Mellin transform does not exist. However, we can solve this problem by rewriting 
the augmented density in a similar way as equation (A.l): defining n as the smallest natural number such that n > fig, we obtain 



I ifr\ p xh(n\( r Y m+2kS ( r 25 



26 \A-« 



For the individual terms the Laplace-Mellin formalism does apply, and the distribution function becomes 



F(G,L) 



Po 
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ds. 



(A. 11) 



(A. 12) 



Now, it is easily observed that we can choose the integration paths C(jq to be identical, as the contours C\, C2 or C3 defined above. 
Therefore the summation can be performed inside the integral. Using the properties 



r(jc + k) = (x) k F(x) and Y(x -*) = (-!)' 



■ Fix) 
(x)k ' 



with (x)k = x(x + 1) ■ • ■ (x + k - 1) the Pochhammer symbol, we find 



s-fio s- f3 a 



+ 1 



n;lj. 



Finally, with the identity 



2 F l (-n,b,c; 1) 



(c - b)„ 

(C)n ' 



(A. 13) 



(A. 14) 



(A. 15) 



the equation for the distribution function also reduces to equation (A. 9). 

If the convergence criterion for path C\ is valid, then every contour can be taken as a line - /oo to Sk + i°° with f3o — k6< 
Sk < Am + in - k)6. Again, for each k the real part of the integrand is positive, so that the distribution function is well-defined. 

Appendix B: A practical series expansion for the distribution function 

We now seek a more practical form for the distribution function. To this aim we first consider the special case where 6 is a rational 
number, denoting 6 = Then we can write 



Fi£,L) = 



Po 



nr(l+p)(8 \ p - V2 1 r r(-> + ns) F(^ M - ns) I L 2 



[ifroj 2m Jc rfl 



M(2^ ) 3/2 n-0s) Wo) 2niJ c r(l-«w)r(p-i+»w) 
Now, using the multiplication formula for the Gamma function 



ds. 



Fikx) = i2nf- k » 2 k kx - 1 ' 2 j~j r(x + Q 



we can write the integral in the form of a Meijer G-function: 

po 1 n n ^-^ )tm F(l 



FiS,L) 



M(27n/r ) 3/2 i2n)"- m mP- 1 ' 2 T( 



-fis) Wo I 



2rl& 



(B.l) 



(B.2) 



(B.3) 



with 



a n+i — 



bn+i — 



n 
i- 1 


H 

+ P - 


m 

-1/2 


m 




m 


i- 1 






n 


m 




i- 1 


m 


- 1 




H 





Since for a general Meijer G-function 



for ; 


= 1,.. 


.,n, 




for ; 


= 1,.. 


■ ,tn, 




for i 


= 1... 


■ ,n, 




for ; 


= 1,.. 


. ,m. 


(B.4) 



(B.5) 
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with p — q and z,a, b real, the convergence criterion for path C\ states that (Mathai 1993) 

1 



m + n- z (p + q) > (B.6) 



or 



m + n--(p + q) = 0, (B.7) 

i p 

2*>i-£fli<-l. (B.8) 

i=0 i=0 

we indeed obtain the conditions 6 < 1 or 6 = 1 and p - | + /3<5 > 0, as used in Appendix A. 

This Meijer G-function can be calculated as a sum of generalized hypergeometric functions (Gradshteyn & Ryzhik 1965, 
eqs. 9.303 and 9.304). For z < 1 the following equation is valid: 



^m+njn+n 



n lYr(*,-* i )fl r ( 1+6 <- fl ') 

: -z 



mbj 



,=1 Y] r(i+bi-bi) Y] r( ai -b t ) 

l=n+l I=n+l 
X m+nFm+n-\ (1 + bj — d\ , . . . , 1 + bj — Cl m+n J 

1 + Z>, , - ft! * 1 + b t - b m+n ; (-l) m+n z m ) , (B.9) 

where the prime by the product symbol denotes the omission of the product when i = I, and the asterisk in the hypergeometric 
function indicates the omission on the i th parameter. Analogously, the equation for z > 1 reads 



Y~[ T(di ~ai)Y\ r (! +*/- «;) 
a \ ;=i 1=1 

| / i m+n m+n 

[[r(i+a,- a,-) |~[ r(a ; - fo/) 



m(a,-l) 



/=«+! 



x m + n F m+n -i (1 + ft; a,-, . . . , 1 + £> m+n a,- ; 

1 + ai -a h . ..,*,..., 1 +a m+ „-ar, (-l) m+ V m ). 
With the coefficients from equation (B.4), we obtain for L 2 < 2r 2 £ 



(B.10) 



m+njn+n 



r2 \ m 



2r 2 <5 



n-l 



n^rr- 

/=o v / /=o 



Poo - A) / + i 

— + 

m n 



m-l 



*z 

7=0 

Now, with the aid of the identities 
and 

we can simplify this expression to 

■2 \ m 



n-l 

n 



1 -j8b + / | i 



m-l 



Poo-Po l + i 
— + 



)m-l / 



7 Z=0 



p+Po-l/2 + l _ i_ 
m n 



p+Po- 1/2 + 1 i 
m n 



n-l 



1=0 



m- 1 



II: , HI 



J 1=0 



m n 



2 \— po+iS+mj 



2rl& 



r(x)T(l -x) 



sin(7rx) 



n 

k=\ 



kn\ n 



n-l 



(B.ll) 



(B.12) 



(B.13) 



* m+n, m+n 



\zr a o/ f / . =0 ^. =Q + ; 



/=() 



+ — +.1 

m n 



i=o \ " ' /=o 

2 \-Pa+i6+mj 



m-l 



l-fi + ^ i 
+ - + 7 

m n 



/=() 



/? + /3 - 1 /2 + Z z 



(B.14) 
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so that, using again equation (B.2), the distribution function reduces to 

\P-3/2 



F(G D Po TQ.+p) le\ p 



n-l oo Y<mp>- + i + jn) I J 2 yPo+iS+mj 

r(i + i + >i»)r(i-jSb + « + /ifi)r(p+jSb-i/2-i<y-7m) v U^s/ 



i=0 7=0 

Finally, the double summation can be grouped into a single index k = i + nj, and we obtain for L 2 < 2r 2 £ 



F(8 L) - Po (*\ V W ( 1| ffl 161 

M(2^ ) 3/2 Uo/ ^\k)T(l-fio + kS)T(p+J3o-l/2-kS)\2riS) ' 

Similarly, for L 2 > 2r 2 £, we find 

Frfi t) po / s f 3/2 y (M r(i +P ) / jM^--" 

M(2^ ) 3/2 Uoj ^U/r(l-/3 M -M)r(p+/3 M -l/2+^)Ur 2 6j 

Although these expressions have been derived for rational values of 5, they can be generalized to any real value, since these functions 
are continuous in 5. Hence we indeed obtain equation (36). 



